## ********************************
##     Figure 3 Subsample Analysis       
##        Hao Zhang & Ye Zhang
##             2025/7/24
## ********************************

# load packages
library(readstata13)
library(lfe)
library(stargazer)
library(tidyverse)
library(doBy)
library(data.table)
library(DescTools)
library(haven)

# set working directory
setwd("../analysis data")
load("firm analysis.RData")

# break down at different ownerships
soe <- subset(firm, soe == 1)
fie <- subset(firm, fie == 1)
private <- subset(firm, private == 1)

# run inverse hyperbolic sine
ihs <- function(x) {
  y <- log(x + sqrt(x ^ 2 + 1))
  return(y)
}

# run regressions
coeff <- function(dataset){
  m1 <- felm(insurance_dummy ~ gdp_pressure + ihs(employment) +
               ihs(wage) + ihs(output1) + ihs(profit) + 
               log(FDI + 1) + gdpgrowth + log(gdp/population) +
               ihs(deficit/gdp) +  ihs(capacity) + begin | 
               frdm + year + educ + ethnicity  + gender + 
               tax + connection + seq | 0 | citycode, data = dataset)
  m2 <- felm(medical_dummy ~ gdp_pressure + ihs(employment) +
               ihs(wage) + ihs(output1) + ihs(profit) +  
               log(FDI + 1) + gdpgrowth + log(gdp/population) +
               ihs(deficit/gdp) +  ihs(capacity) + begin | 
               frdm + year + educ + ethnicity  + gender + 
               tax + connection + seq | 0 | citycode, data = dataset)
  coeff <- matrix(NA, 2, 3)
  coeff[1,] <- c(coefficients(m1)[1], confint(m1, level = 0.95)[1,1:2])
  coeff[2,] <- c(coefficients(m2)[1], confint(m2, level = 0.95)[1,1:2])
  return(coeff)
}

result <- matrix(NA, 6, 3)
result[1:2,] <- coeff(soe)
result[3:4,] <- coeff(fie)
result[5:6,] <- coeff(private)

# draw the final graph
pdf("gdp_pressure_ownership_firm_dummy.pdf")
par(mfrow = c(1,1))
par(mar=c(5,6,3,2))

# plot ownership
plot(NA, xlim = c(-2, 0.5), ylim = c(1, 11), yaxt = "n", 
     xlab = "Effect of GDP Growth Pressure", ylab = "", cex.lab=1.5, cex.axis=1.2)

axis(2, at=3, labels="State", las=1, cex.axis = 1.3)
axis(2, at=6, labels="Foreign", las=1, cex.axis = 1.3)
axis(2, at=9, labels="Private", las=1, cex.axis = 1.3)

for (i in 1:3) {
  points(y= 3*i + 0.5, x= result[2*i-1, 1], pch = 16, cex = 1.4, col = "grey30")
  points(y= 3*i - 0.5, x= result[2*i, 1], pch = 17, cex = 1.4, col = "grey30")
  segments(y0 = 3*i + 0.5, x0 = result[2*i-1, 2], y1 = 3*i + 0.5, x1 =  result[2*i-1, 3], lwd = 1.8, col = "grey30")
  segments(y0 = 3*i - 0.5, x0 = result[2*i, 2], y1 = 3*i - 0.5, x1 =  result[2*i, 3], lwd = 1.8, col = "grey30")
}

legend("bottomleft", legend = c("Unemployment", "Medical and Pensions"), col = c("grey30", "grey30"), pch = c(16, 17), cex = 1.2)
abline(v = 0, col = "Red")

dev.off()
